Analysis of the discontinuous Galerkin method for elliptic problems on 

surfaces 
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We extend the discontinuous Galerkin (DG) framework to a linear second-order elliptic problem on a 
compact smooth connected and oriented surface in R . An interior penalty (IP) method is introduced on 
a discrete surface and we derive a-priori error estimates by relating the latter to the original surface via 
the lift introduced in |Dziuk|jl988^ . The estimates suggest that the geometric error terms arising from 
the surface discretisation do not affect the overall convergence rate of the IP method when using linear 
ansatz functions. This is then verified numerically for a number of test problems. An intricate issue is the 
approximation of the surface conormal required in the IP formulation, choices of which are investigated 
numerically. Furthermore, we present a generic implementation of test problems on surfaces. 

Keywords: discontinuous galerkin; interior penalty; surface partial differential equations; error analysis. 



1. Introduction 

Partial differential equations (PDEs) on manifolds have become an active area of research in recent years 
due to the fact that, in many applications, models have to be formulated not on a flat Euclidean domain 
but on a curved surface. For example, fliey arise naturally in fluid dynamics (e.g. surface active agents 
on the interface between two fluids, 'James & Lowengrub'('2004')) and material science (e.g. diffusion of 
species along grain boundaries, Deckelnick et al. (2001)) but have also emerged in areas as diverse as 
image processing and cell biology (e.g. cell motility involving processes on the cell membrane, |Neilson] 
etal. (201 l|l or phase separation on biomembranes, 'Elliott & Stinner (2010l). 

Finite element methods (FEM) for elliptic problems and their error analysis have been successfully 



applied to problems on surfaces via the intrinsic approach in Dziuk ( 1988 1 based on interpolating the 
surface by a triangulated one. This approach has subsequently been extended to parabolic problems in 
Dziuk & EUiottl ( |2007a| l as well as evolving surfaces in [Dziuk & Elliott| ( [2007b| . |Ju & Du| ( p009l l and 
Giesselmann & Mtiller| (|2013 ) have also considered finite volume methods on surfaces via the intrinsic 
approach. However, as in the planar case there are a number of situations where FEM may not be the 
appropriate numerical method, for instance, advection dominated problems which lead to steep gradients 
or even discontinuities in the solution. 

DG methods are a class of numerical methods that have been successfully applied to hyperbolic, 
elliptic and parabolic PDEs arising from a wide range of applications. Some of its main advantages 
compared to 'standard' finite element methods include the ability of capturing discontinuities as arising 
in advection dominated problems, and less restriction on grid structure and refinement as well as on the 
choice of basis functions. The main idea of DG methods is not to require continuity of the solution 
between elements. Instead, inter-element behaviour has to be prescribed carefully in such a way that 
the resulting scheme has adequate consistency, stability and accuracy properties. A short introduction 



to DG methods for both ODEs and PDEs is given in Cockbum (2003 1. A history of the development 
of DG methods can be found in |Cockburn et a/.| ( |2000| l and [Arnold ef a/.| ( |2002| l. [Arnold ef oT] ( |2002[ l 
provides an in-depth analysis of a large class of discontinuous Galerkin methods for second-order elliptic 
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problems. 

The motivation of this study has been to investigate the issues arising when attempting to apply 
DG methods to problems on surfaces. We restrict our analysis to a linear second-order elliptic PDE 
on a compact smooth connected and oriented surface. We expect that parabolic problems on evolving 



surfaces as featured in the above mentioned applications can be dealt with along the lines of Dziuk & 
|EUiott| ( |2007b| l. 

This paper is organised in the following way. We consider a second-order elliptic equation on a 
compact smooth connected and oriented surface F C M^ and introduce a particular DG method known 
as the interior penalty (IP) method on a triangulated surface f/,. The surface IP method we consider 
is similar in nature to the one introdu ced in|Arnold|(|1982[ ), and its well-posedness follows naturally 
from results in the planar case given in Arnold et al. ( 2002[ l and Ainsworth & Rankin ( 2011[ l. We then 
derive a-priori error estimates in the appropriate norms by relating f/, to F via a lifting operator and by 



making use of results from Dziuk ( 1988 1 and Giesselmann & Miiller (20131) to show that the additional 
geometric error terms arising when approximating the surface scale in such a way that they do not affect 
the convergence rates proved and observed for the standard FEM approach in |Dziuk| ( |1988| l when using 
linear ansatz functions. 

We then present some numerical results, making use of the Distributed and Unified Numerics Envi- 
ronment (DUNE) software package (see Bastian et al. (2008b|, Bastian et al. ( 2008a I) and, in particular, 
the DUNE-FEM module described in Dedner et al.^ i ^lQlQ) (also see dune.mathematik.uni-freiburg.de 
for more details on this module). We consider a number of test problems, for which we compute ex- 
perimental orders of convergence (EOCs) in both the !?■ norm and the DG norm, and show that these 
coincide with the theoretical error estimates derived in the previous section. Furthermore, we consider 
several intuitive ways of approximating the surface conormal in our IP formulation, and investigate the 
resulting schemes numerically. In the process, we present a generic implementation of test problems on 
surfaces which follows as a du^ect application of the |Demlow & Dzri3c, ( |2008) algorithms. 

Finally, we briefly present numerical results for nonconforming grids and higher order polynomial 
ansatz functions, which suggest that the convergence rates of the standard FEM approach still hold for 
such generalisations. 



2. Notation and Setting 

The notation in this section closely follows the one used in |Dziuk| ( |1988| l. Let F be a compact smooth 
connected and oriented surface in M? . For simplicity, we assume that dF = 0. Let d denote the signed 
distance function to F which we assume to be well-defined in a sufficiently thin open tube U around F . 
The orientation of F is set by taking the normal v of F to be in the direction of increasing d whence 

v(<^) = Vrf(^), ^GF. 

With a slight abuse of notation we also denote the projection to F by ^, i.e. ^ : f/ — > F is given by 

^{x)^x — d{x)v{x) where v(;c) := v(^(.:c)). (2.1) 

Later on, we will consider a triangulated surface Fi, cU approximating F such that there is a one-to-one 



relation between points x e f/, and ^ e F so that, in particular, the above relation ( 2. 1 1 can be inverted. 
Throughout this paper, we denote by 



P(^):^/-v(^)®v(^), ^GF, 
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the projection onto the tangent space TpF on F at a point £, E P. Here (g) denotes the usual tensor 
product. 

Definition 2. 1 For any function tj defined on an open subset of U containing F we can define its 
tangential gradient on F by 

Vri? := V77 - (Vt] • v) V = PVrj 

and then the Laplace-Beltrami operator on F by 

4rrj:=Vr-(VrTj). 
Definition 2.2 We define the surface Sobolev spaces 

H"\F) := {u e L^{F) : D"u e L^{F) V|a| s^ m}, m e NU{0}, 
with corresponding Sobolev seminorm and norm respectively given by 

\\a\=m / \k=0 / 

We refer to |Wloka] ( |1987| l for a proper discussion of Sobolev spaces on manifolds. 

The problem that we consider in this paper is the following second-order elliptic equation: 

-Aru + u = f (2.2) 

for a given / e L^{F). Using integration by parts on surfaces the weak problem reads: 

(Pr) Find u e H^ (F) such that 



Vru-Vrv + uvdA= / fvdA Vve//'(F). (2.3) 

r Jr 

Existence and uniqueness of a solution u follows from standard arguments. We assume that u e H^{F) 
satisfies 

MHHr)^C\\f\\^2^r) (2.4) 

where we refer to Aubin ( 1982|l and Wloka ( 1987[) for more details on elliptic regularity on surfaces. 



3. Approximation and Properties 

To obtain a discretisation of u, the smooth surface F is approximated by a polyhedral surface Ff, C U 
composed of planar triangles. Let 5^, be the associated regular, conforming triangulation of iT, i.e. 

J7,= U ^h- 

The vertices are taken to sit on F so that f/, is its linear interpolation. We assume that the projection map 



^ defined in (2.1 1 is a bijection when restricted to I/,, thus avoiding multiple coverings of F by I/,. Let 



?/, denote the set of all codimension one intersections of elements K^ ^K^^ G .% (i.e., the edges). We 
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define the conormal n^ on such an intersection e/, e d}, of elements Kj^ and Kj^ by demanding that 

• n^ is a unit vector, 

• n,| is tangential to (the planar triangle) Kj^, 

• in each point x e e/, we have that n'^ ■ (y — x) ^0 for all y e /T^^ . 

Analogously one can define the conormal n^ on e/, by exchanging K^ with /T^J" . A discrete DG space 
associated with F/, is given by 

V,, := {vv, e l2(I7,) : v„|^,_ G P'(/^,,) V/T/, e ^-J 
i.e. the space of piecewise linear functions which are globally in L^(I/,). For v/, e V/,, let 

"h -''h\dK:'-- 

We can now define a discrete DG formulation on f/, for a given function //, e L^iPh) (note that, in gen- 
eral, this is not a finite element function, it will be related to the function / given in problem (Pj-) later 
on, see p.5|l below): 



(Prj) Find m/, £ V/, such that 

a^{uh,Vh)= Y^ / fhVhdAh'^VheVh (3.1) 



^;,G^, ^'' 



where 

JPi 



K,,e.%-'^h 



a'£{uh,Vh):^ Y, / '^r,Uh-'^r,Vh + UhVhdA 



h 



/■ 1 1 

+ L I ^eM~H)i^t-^l)dsh- (3.2) 

The penalty parameters jSe,^ are given by jS^,^ = (Oefh'^^ where he^^ is some length scale associated with 
the intersection e/, (for instance, the edge length). The interior penalty parameters (Oe,^ are uniformly 
bounded with respect to /i := max^^^g^^^ /ig,^. 



Remark 3.1 This formulation corresponds to the one found in Arnold et al. (2002i in the case when 
the do main is flat and is similar in nature to the original formulation of the IP method found in Arnold 
( 1982 1 for which the conormals n^^ ' are associated with their respective gradient terms rather than the 



scalar terms. It is important to point out that this formulation is not equivalent to using the formulation 
found in [Arnold et al.\^002\ on f/,. We will discuss this issue further in Section[5] 

We now define a norm on the space of piecewise smooth functions: 

Definition 3.1 For m/, e V/, we define 
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The DG norm is given by 

Lemma 3.1 Let c^k^ denote the set containing the individual edges of element A),. Then if jSe,^ = C0e,,hei} 
with 

1 V- leAp 
COe^> max - 2^ — — forall e/, e^/,, (3.3) 

ei,cdKi, '"^*'^A 

then Up" is stable and bounded. Hence there is a unique solution m/, e V/, of (Pj^) which satisfies 

\\uh\\DG^C\\fi,\\^2(r„)- (3.4) 



Proof. Boundedness and stability of (3.2i follow in a similar way as for the classical IP method (see 



Arnold et al. ( 2002|l for more details) since all the arguments apply to f/,. For the lower bound of the 



penalty parameters, the proof of Lemma 2.1 in |Ainsworth & Rankin | ( [201 1[ ) applies straightforwardly 

to the surface f/,. Note that the reason why these results naturally extend onto f/, is that the latter is 

composed of planar triangles. By Lax-Milgram, the uniqueness property follows. D 

Our goal now is to compare the solution u e H^{r) of (Pj-) with the solution Uh G V), of (Pp") but 



these functions are defined on different domains. The approach suggested in Dziuk (1988il is to lift 
functions defined on the discrete surface f/, onto the smooth surface F. 

Definition 3.2 For any function w defined on f/, we define the lift onto F by 



where by ( 2. 1 1 and the non-overlapping of the triangular elements, x(^ ) is defined as the unique solution 
of 

Extending w' constantly along the lines s t-^ ^ +sv{^) we obtain a function defined on U. In 
particular, we 

define //, such that /,( = / on F. (3.5) 



By (2.1 1, for every Ki, G =5/,, there is a unique curved triangle KJ^ :— t, {K/,) C F. Note that we assumed 
£, (x) is a bijection so multiple coverings are not permitted. We now define the regular, conforming 
triangulation 3^/ of F such that 

r= U <■ 



<e.-^,' 



The triangulation 5^^' of F is thus induced by the triangulation 5^, of I/, via the lift. Similarly, ej^ ;= 
^ (e/,) € (o/^ are the unique curved edges. 

The appropriate function space for lifted functions is given by 

Vl := {vi e L^{F) : yj,(^) = v,,{x{^)) with some y,, G V,,}. 



Note that the DG norm for functions wj^ G V^ is the same one as in Definition 



3.1 



but with the triangula- 
tion £^,' instead and corresponding length scale hi associated with e', . The context of its use makes it 

h 

clear which DG norm we are dealing with. Furthermore, we observe that 

h'e,, ^ he,, (3.6) 
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since the deformation of the straight edges can only increase their length. We now prove some geometric 
error estimates relating T to i^. 

Lemma 3.2 Let F be a compact smooth connected and oriented surface in M? and f/, its linear inter- 
polation with outward unit normal V/,. Let H — V^d and Pi, — I — V/, (^ Vh- Furthermore, we denote 
by 5/, the local area deformation when transforming K^ to /T^ i.e. 5hdAf, = dA and Se/^ the local edge 
deformation when transforming e/, to ej^ i.e. Sg^dsf, = ds. Then we have 

\\d\\L^(r)^Ch^, ||l-5,,||i»(r) ^C/!^, ||v- v,,||z,^(r) s^ O, \\P -Rh\\L~(r) ^Ch^ 
\\l-de„\\L-{r)^Ch\ ||P-/?,J|i.(^)S=:C/z2and||«-P«j,||^.(^,)<C/22 

where/?/, := ^P{I - dH)Ph{I - dH) and/?^,, := ^P{I -dH)Ph{I-dH). 

Proof. All these geometric estimates follow from standard interpolation theory via the linear interpola- 
tion of F. Proofs of the first four estimates can be found in Dziuk ( ,1988) , the fifth (and thus sixth) one 
in Giesselmann & Miiller (2013 1. The last estimate is a corollary of another result in Giesselmann & 



Miiller 



(12013'), which states that if n and n'l^ are given as before and T denotes a unit tangent vector on 



some ejj e (f^, we have 



\{T,n[)\^Ch\\\~ {n,n[)\^Ch^ 



Writing fn'^ = (T,f njjT+ {n,Pn'i^)n, we deduce that indeed 



\n-Pn' 



■hh-iel) 



= \\n-iT,Pn',)T^{n,Pn[)n\ 

< |1 - {n,Pni)\ + \i^M)\ = |1 - («,4)l + 1(^,4)1 = 0{h'). 



L-{e') 



Lemma 3.3 Let m/, e Vh satisfy (3.4 1. Then u[ e V/ satisfies 



for sufficiently small h. 

Proof. We first show that that for any function v/, G V/,, 

||v/,||dg^C||v{||dg- 
The I • I J ;j component of the DG norm is dealt with in exactly the same way as in 



D 



(3.7) 



(3.8) 



3.2 



and the fact that h^ '^ h , (see (3.6l), we obtain the following for the 



Dziuk 



( 1988 1. Similarly, 



making use of Lemma 
component of the DG norm 

Lk^I H~-,fdsh-;^ ILK' Is 



\*ji 



h '^^h 



'+-v'-? 



1 ds 
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which yields the desired estimate for sufficiently small h. Noting that ||//,||£2(j-) ^ ^WflWi^ir) = 
C||/||£2(j-) (see Dziuk ( 1988 1), we can extend the stability estimate (3.4 1 to the lifted discrete function 
wjj as required. 



D 



We now define a bilinear form on F induced by a^^ and the lifting operator: 



Op (m,v) 



f \ 1 

e' eS' "a ^ ^ 



£ ( Bj{u+-u-){v+-v-)ds 



(3.9) 



e'eS' *< 



w+ 



where n+ and n are respectively the unit surface conormals to K'l^ and A",' on ej^ € ^^' and the penalty 
parameters are defined to be j3 ; '-=4^- This bilinear form is well defined for functions m, v e H^{r) + 
V}', and since the weak solution u given by ( |2.3| ) is in H^{r) it satisfies 



(3.10) 



«r («,v') = £ / /v rfA Vv e H\r)+Vl 



We now state and prove a technical estimate of importance for boundedness and stability of Jf . 
Lemma 3.4 Let w e //^(r) and wj, e ¥";(. Then for sufficiently small h. 



,'m|2 



l!Vr(w + wl)||^,(^^,, ^C -|lVr(w + vvi)|l^,(^,j+/z|lw||^,(^,, 



(3.11) 



Proof. We define vv and w/j such that their lifts coincide with w and w'j^, respectively. Since w + Wh £ 
H^{Kij) on each Ki,, applying the trace theorem and a standard scaling argument on K^ e J^ yields 

JdK,, \hJKi, Jk,, '' 

where we used that V^ w/, = thanks to the linearity of the finite element functions. Lifting the estimate 
onto r and making use of estimate (2.17) in |Demlow| ( |2009[ ), we have 



< 



Vr{w + wi) -Reyriw + wi) ds ^C ( y f \Vr{w + yv'i^)\^ dA + h f \Vlw\^ + \V rw\^ dA] 

V« J 4 Jk' J 



where Re^ is given as in Lemma 3.2 We thus obtain 



(l-Ch^) [ \Vriw + w[)\^dss^c(j[ l^^riw + wi)]^ dA+h f \Vlw\^ + \Vrw\^ dA 
which yields the desired inequality for h small enough. 



D 
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Following the lines of Arnold et al. (20021 along with Lemma 3.4 we can show the boundedness 
estimate 

a'f{w + wl^i)^C',{\\^v + ^v[\\DG + h^\MHH^))h'k\\DG forall w^H\r), w^v^^vi (3.12) 

and the stability estimate 



«rK,w;j^Ql|wi||^G forall w^ V,' 



(3.13) 



4. Convergence 

Theorem 4.1 Let u e H^{r) and m/, e V/, denote the solutions to (Pj-) and (Pj^), respectively. Denote 
by M/j G y^J the lift of m/, onto F. Then 

ll«-«/,llL2(r)+/^ll«-4lloGS^C/!^||/|lL2(r)- 



The proof will follow an argument similar to the one outlined in [Arnold et al. ( 2002| l. Using the 
stability result (|3.13 i, we have 



C^. L,j L-j 

where 0^' G V}'. Since we do not directly have Galerkin orthogonality the first term is not zero, and the 
second term will require an interpolation estimate. These terms are addressed by the following lemmas: 

Lemma 4.1 For a given w e H^{r) there exists an interpolant I^w G V^ such that 

\\^-t'hM\LHn +''\\^r{w-llw)\\L2^r) ^ Ch^ (l|Vf vHlz,2(r) +/^l|Vrvv||i2(r)) ■ 
Proof. See |Dziuk| ( |T988l ). D 



Lemma 4.2 Let u and mJ^ be given as in Theorem 4. 1 and define the functional £/, on V^ by 
Then £/, can be written as 



L /,K+-«/r)^(Vrv{+-«+-Vrvi--n 



il«W 



l-^)K^A 






,/+ .J4 



,.;- .J--\ 



« - «;r ) 2 (^/T (^ - dH)PVrv'+ ■ n'+ - P, (/ - dH)PVrv',- -n',-)^ 



ds 



I /,K+-v{r)^(Vr4+-«+-vr«jr-«-) 






./- „/-^ 



(^r - ^r ) 2 « ('f - dH)PVru',: ■ «;+ - P,7 (/ - dH)PVru',; -n, ij 



ds 
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where /?/, is given as in Lemma [3T2| Furthermore, £/, scales quadratic ally in h i.e. 

mv'„)\^Ch^f\\^2,nH\\DG- 



(4.2) 



Remark 4.1 Note that the error functional £/, in Lemma \4^ includes all of the terms of the classical 
FEM setting (see Dziuk| ( 1988| l) as well as additional terms arising from the jumps across elements 
which characterise the DG method. 



The proof of Lemma 4.2 will be the main part of this section. Before we give its full proof, we will 



complete that of Theorem 4.1 assuming this result. Using the estimate (4.1 1 given at the start of the 
proof of Theorem |4.1| the boundedness result ( |3.12| i, the elliptic regularity result p.4| l and the quadratic 
scaling of £/, in h (4.2 1, we have 



Uh - u'hWlc < ^Ef,{(t)i - mI) + ^a'f{(t>l, ~u,<pl~ 4) 



€ 



\DG 



thus 



-iEh{<pl-u[) + -^^(\\^i-u\\DG+h^\\u\\HHr))\Wh-ui\\i 

<C/l2||/||^2(^)||,^^-«;,||oG+c(||(/)i-M||oG + /l'll/llz,2(^))||0i-«i||DG, 

Ul-u[\\DG^Ch^\\f\\LHr)+C\Wh-u\\DG- 



4.1 



we obtain 



Now taking the continuous interpolant 0^' ~ I^u and using Lemma 

\\u-4\\DG^\\u-^j,\\DG + Ui-ui\\DG^\\u-^h\\DG + Ch\f\\^2^r)+C\\K-4^^ 

as required. The L^ error estimate can be derived using the usual Aubin-Nitsche trick in a similar way 
as in Dziuk (1988|l, which concludes the proof of Theorem 4.1 



Proof of Lemma \4^ The expression for the error functional £/, given in Lemma K2l is obtained by 



considering the difference between the two equations (3.10i and (3.1 1. In order to do this, the integrals 
of p.l| i have to first be hfted onto F. For every K/, e =3^,, we have 

/ Vj-ftM/! • V/7, v/, + M/,y/, dAh = / Rh^ru[ ■ ^rA + -^u'hv[ dA. 
JKi, Jk[ Oh 

Furthermore, for every e/, e (I/,, we have 

r I 1 



,./- 



+ i^k~0^iKiI-dH)P'7ru 



"h >2^'h 
And finally, we have using j3^/ = -^ that 



'+ J+ 



Pf;{I-dH)PVrii 



I- „/-^ 



ds. 



/ h, (4 - «/, ) i^'h - ^'h ) dsi, = ,pj («^ - "1 ) (^'^ - ^i ) ^^• 

Je^ Je> I' 



lOofEB 



The right-hand side of ( 3. 1 1 gets transformed in a similar way 

fhVh dAi, = 






K>G.'7'-"'l 



fv[^dA. 

Oh 



Making use of the above, the difference between the two equations p.lO[ i and (3.1 1 yields 
= flf(M,vj,)- ^ lfv[dA-a'f^^{uh,Vh) 









= a'r^{u-ui,vi)^Eliv',) 



as required. 

Finally we need to show that the error functional Ef, scales quadratically in h i.e. 

mvi)\^Ch\f\\^2^r)hi\\DG- 
To this end we need to show that the additional terms arising in the error functional £/, do not affect 



the convergence rates expressed in Dziuk ( 1988 1. The first term of the error functional £/, (the element 



integral) is the one resulting from the standard surface FEM approach. By Lemma 3.2 this term scales 



quadratically in h and making use of the stability estimate (3.7 1 this term scales like the right-hand side 
of (|4.2|l. We will now get a bound for the third term of £/,, for which we have the following: 



eiG4 



/+ 
.1^ h ' 



i< 



2(Vr« 



^•«"- 



r«l •« )(l + ^ 



-VrU- 



^ei, 



1 



ff,, 



i-: 



/,+ " ^'/r ) I (Pk it - dH)PVru'+ ■ n'+ - P, (/ - dH)PVru\; ■ n\; ) ^ ds 



' _ „i Je' ^ \ Oe 



4e<"* 



1 



(v-r- 



,./- 



)^((Vr<+-«^ 



VrW 



ruh 



iP+{I-dH)PVru 



1+ J+ 



Pi;{I~dH)PVru\ 



ds. 



Making use of standard arguments as found in [Arnold et a/.| ( p302[ l along with Lemma 3.4 Lemma 3.2 
and the stability estimate (3.7i it is clear that the first component in the above scales appropriately, so 
all we have to deal with is the second component. We first note that 



VrMj,+ -«^ 



P+{I~dH)PVru 



VrW 



\:-n^ 



'Vru'+-P{I-dH)n'+ 



/+.,/+ 
<' 



Vru'+-n+- 
yru'+ ■ («+ 



Vru'+-P{I-dH)P+n 



Pn 



l+\ 



-dHVru'J-n' 






h ' 



hence 



e'8. 



M 



-)\{{^ru[+-n+-Vru\;-n-) 



{P+{I~dH)PVru\ 



1+ J+ 



-P^{I-dH)PVru'i; 



e'5e 



-(-: 



-)\{in^ 



ds 



-Pn'+) ■ Vru'++dH{Vru[+-n[+ - Vpu',; ■ n'^) 



\P< 



Vj-M^ I ds 
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For the first component of the above, we have 



,; ^gl •I",, <^e,, ^ 



^4e< '' ^ "■ 



after applying Cauchy-Schwartz. Using similar arguments as done for proving boundedness of the 
classical IP method (see [Arnold et a/.| ( |2002| l), we have 






.+ D,J+ll2 11.,/ 112 



<Cmax|l«+-Pn|,+ |l2 £ ^ /; , || V^m' H^ ,^ Cmax |1«+ -PwJ+n;^, , j|m}. 



'h^H " K[e.3rie\^edKl '^" " 'h'^H 



esl ^ ^'A^ 



where the last inequality was derived using similar arguments as in the proof of Lemma |3.4| For the 
second component, we have 



I / f <■ 



{^\:-^];)\{dHVru]^-n[+)ds 



1 1 . /...„ ,+ ,+\2 



1 



Pursuing the analysis as before and using the fact that the Hessian H is symmetric and bounded, we 
have 






4e-^^;,e^4 



where again the last inequaUty follows from applying similar arguments as in the proof of Lemma 3.4 
We can now estimate the error functional £/,: 

+ l|l-^llL"(r)ll/llL2(nlkAllDG + Cmax||«-P«{,||^^,,J|M[llDGl|v'/J|DG 

+ C||«'||L-(r)||4l|z3G||vi||DG. 
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So by Lemma 3.2 and the stability estimate ( |3.7[ i we have 

\Eh{v',:)\^Ch\f\y^r^\\vi\\oG 
for every v'^ e V^ as required. 



n 



5. Numerical Tests 

5.1 Implementation Aspects 

The IP method has been implemented using DUNE-FEM, a discretization module based on the Dis- 
tributed and Unified Numerics Environment (DUNE), (further information about DUNE can be found 
in |Bastian et a/.|p008al ), |Bastian et a/.| ( |2008b| l and |Bastian et a/.|P012l i). DG methods are well tested 
for the DUNE-FEM module, as shown in 'Dedner et al. (20101, Brdar et al. (12012), but only simple 



schemes have been tested for surface PDEs (further information about the DUNE-FEM module can 
be found in |Dedner et al.\^20\0\ and [Dedner et al.\^20\2\ ). \n all our numerical tests we choose the 
polynomial order on each element K^ € 5^, to be 1 and interior penalty parameters to satisfy ( 3.3 1. The 
initial mesh generation for each test case is performed using the 3D surface mesh generation module of 



the Computational Geometry Algorithms Library (CGAL) (see Rineau & Yvinec (2009 1). 

When performing mesh refinements it is often the case that there is no explicit projection map 
for mapping newly created nodes from f/, to F, hence t,{x) must be approximated. Two different 
algorithms, discussed in more detail in Demlow & Dziuk| ( |2008 1, have been tested for such problems: 
one being Newton's method and the other being an ad-hoc first-order method. Assume that xq E U and 
that we wish to compute ^{xq). The Newton method seeks to find a stationary point of the function 
F{x,X) = \x — xq\^ + X(J){x) with starting values chosen to be (xo,Ao) ~ (xo,2(^(xo)/|V(^(xo)P), where 
is the level-set function of F (and not necessarily a signed-distance function). We iterate the method 
until 



0W' 



V0(x) 



-XQ 



|v</)W| 



-xo\ 



1/2 



<tol 



(5.1) 



is reached. Note that this stopping criteria incorporates both how close the iterate is to the surface F 
as well as how accurately it lies in the normal direction. The second method is given by the following 
first-order algorithm: 

1. Stipulate tol and xq and initialise x = xq- 

2. While ( |5.1| l is not satisfied, iterate the following steps: 

(a) Calculate! — x— lyL 1 2 and dist — sign{^ {xo))\x ~ xo\. 

(b) Set x^xo~ dist ^^^ ■ 

We can make this algorithm more flexible by only requiring a second order finite difference approx- 
imation of V(j){x). It was observed in Demlow & Dziuk (20081 that in practice the second of the two 
algorithms was more efficient due to the fact that each step of Newton's method is relatively expensive. 
This was observed in all of our numerical tests. It has also been noted that there has not been any rig- 



orous error analysis done for either of the two algorithms with respect to the stopping criterion (5.1 1. 
Though we do not provide any such error analysis, our numerical tests suggest that both algorithms may 
stagnate at certain points for which the stopping criterion tolerance is never reached, even with fairly 
refined initial meshes. The normal direction error contribution of the stopping criteria appears to be 
responsible for the algorithm stagnating, hence for such points we remove this contribution so that the 
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algorithm terminates and the resulting point lies on the surface nevertheless (albeit not necessarily at 

^(xo)). 

In addition, we make use of this algorithm to provide a generic implementation of test problems on 
surfaces. Computing the Laplace-Beltrami operator of some given function over an arbitrary compact 
smooth connected and oriented surface given by the zero level-set of some function is tedious and 
requires changing the implementation for every such surface. In particular, we would need to explicitly 
compute the outward unit normal of the surface and its gradient whenever we consider a new surface. 
For any u £ C^{M?), we have 



Aru^Au~v-V^uv~tr{Vv)Vu-v (5.2) 

where A is the usual Euclidean Laplace operator in M?, V^u e M?^^ the (Euclidean) Hessian of u, Vu 
the (Euclidean) gradient of u and finally tr(Vv) the trace of Vv where Vv G M?^^ whose entries are 
the (Euclidean) partial derivatives of each component of the normal. We can make use of the ad-hoc 
first-order algorithm described previously to approximate the outward unit normal v of F in ( |5.2| l: this 
is done by computing v{£,{xq)) w sign{(p{x())){£,{x()) ~xq) where £,{xq) is the approximation of £,{xo) 
resulting from the algorithm . We may also approximate the (diagonal) entries of Vv via second-order 
finite difference approximations as done for the approximation of V0 in the first-order algorithm. Such 
a generic implementation has the benefit of only requiring input of the level-set function for the surface 
and nothing more, significantly facilitating numerical tests. The error caused by our approximation of 
the Laplace-Beltrami operator appears not to affect the resulting convergence order for any of our test 
cases. 

5.2 Approximation of Surface Conormals 
Consider the IP bilinear form d'f given by 

/■ 1 1 

where n+ and n^ are simply vectors which lie on the intersection e/, G d}, of neighbouring elements K^ 
and Kf^ . Now assume that we want to assemble the system matrix on an element K/j and we assume that 
Ki, = K^ for all e/, C dKi,. To this end, we fix v/, = (p^ with supp(^^) = Ki, which leads to 



JKi, 



/■ 1 1 



E / ^e,XK^"h)^ dSh- 



eHCdKi,-""!' 
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To assemble the block on the diagonal of the matrix we need to take m/, = \j/^ with supp(i//^) = K/,. 
For the off-diagonal block we take m/, = 1//+ with supp(v/+) = Kj^ for one neighbour K^ of A/,. We 

will then discuss different choices for ng/ which are Unked to several intuitive ways of approximating 
respectively the surface conormals n+' . We use one choice for n+ in both cases. To cover all of the 
choices we want to consider, it is necessary to use different choices for n^ , i.e., the vector belonging to 
the element Ki, on which we are assembling the matrix. For the diagonal block we will denote our choice 
for this vector with n^ and use the original notation n^ for the choice used to assemble the off-diagonal 
block. Note that n^ = nj^ for all of the choices discussed below except for Choice 3. 
Now consider m/, = )/^ with supp(v/^) = Kij in ( |5.3| l using n^ instead of n^ : 



JKi, 



f \ 1 

Next we take m/, — i/+ with supp ()//+) ~ Kj^ for one neighbour Kj^ of Kf,, we now have 

/■ 1 1 

We can now prescribe choices for the vectors «^, nj , «+ and will later investigate the behaviour of the 
numerical scheme ( |5.3[ ) for different choices of these three vectors. 
Choice 1 

Such a choice corresponds to using the IP method in a planar setting, for which ny| = —nj^, and is the 
simplest scheme to implement. 
Choice 2 



This choice yields the numerical scheme (3.2 1 that has been discussed up to now and used in the error 
analysis. 
Choice 3 



: ne,. = TYT^ TVT : « 



" 15^-4)1 '• l5K-4)l " \ki<-n,)\ 

This choice corresponds to prescribing the vectors to be the average of the two conormals and yields 
additional symmetry in the resulting matrix due to the fact that the vectors are now independent of the 
element on which they are computed. 
Choice 4 

«D = «/7 > «7, = -«/! > <, = -«A- 

This particular choice corresponds to using the formulation of the IP method found for example in Arnolc^ 



: 
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et al. ( 2002[l on I/,, but with a modified penalty term that does not depend on the conormals i.e. 



a"' 



%{uh,vh) = ^ / Vr^M/, • yr,vh + UhVh dAh 

f I 1 

+ Y^ Pe^ {lit ~ "/7 ) i^ ^ ^'h ) ^^h (modified penalty temi). 
We summarise the choices in Table [T] 



Choice 


«D 


««„ 


«i 


Description 


1 
2 

3 

4 


I^K-«,T)I 


"a 

IH«;-«,T)I 


-"a 

"/I 

IW-«,7) 

ii(«,r-«,7)i 

-"A 


Planar (non-sym) 
Analysis (sym pos-def) 

Average (sym pos-def) 


1 Arnold et al.\\2iX)2) (sym pos-def) 



Table 1: Choices of n^, n+ and n^ , description of the numerical schemes they respectively lead to 
and properties of resulting matrix. 



We also consider the Arnold et al. ( 2002 1 formulation with its true penalty term given by 



E / ^<^h i'-'h'^h +"h'^h)- (^/l"/T + ^h "a ) '^'^h (true penalty term). 

e„CdK,/^h 

Choosing v/j = (p^ and m/, = y/^ as before yields 

ei.CdKi/'^l- 



For M/, = V/+ we now have, 






The matrices arising from Choices 2-4 are symmetric positive definite, so the Conjugate Gradient 
(CG) method is particularly well suited for such matrix problems. Choice 1 however yields a non- 
symmetric matrix, for which we use the Biconjugate Gradient Stabilized (BICGSTAB) method. All 
of these solvers make use of the algebraic multigrid algorithm (AMG) preconditioner coupled with the 
incomplete-LU factorisation preconditioner to speed up the solvers. Information on the implementation 
of these solvers and preconditioners in DUNE can be found in Blatt & Bastian ( 2007| l and on their 
parallelisation in |Blatt & Bastian| p008 ). 

We first tested our code on a sphere where the projection algorithm described in Section [5T| is not 
required. The results showed that the expected convergence rates and the choices of n^, n^ and «+ had 
little influence on the results. Hence we have decided not to include these tests since no insight can be 
gained. 
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5.3 Test Problem on Dziuk Surface 



The first test problem, taken from Dziuk ( 1988 1, considers 



' Apu - 



/ 



(5.4) 



on the surface F = {x e 



(xi — X3) +^2 +X3 = 1} whose exact solution is chosen to be given 



by u{x) — x\X2. The outward unit normal to this surface is given by v(jc) — (x\ —x'^,X2,xt,{1 — 2(jci — 
■^3 )))/(! +4x3(1 —xi — Xj))''^. There is no explicit projection map for mapping newly created nodes 



to F so (^ (x) has to be approximated via the ad-hoc first order algorithm described in Section 5.1 



Elements 


h 


L2 -error 


L2-eoc 


DG-error 


DG-eoc 


92 


0.704521 


0.243493 




0.894504 




368 


0.353599 


0.0842372 


1.53 


0.490805 


0.87 


1472 


0.176993 


0.0268596 


1.65 


0.263808 


0.90 


5888 


0.0885231 


0.00637826 


2.07 


0.135162 


0.97 


23552 


0.0442651 


0.00171047 


1.90 


0.0685366 


0.98 


94208 


0.022133 


0.000416366 


2.04 


0.0343677 


1.00 


376832 


0.0110666 


0.000104274 


2.00 


0.0171891 


1.00 


1507328 


0.0055333 


2.60734e-05 


2.00 


0.0085935 


1.00 



Table 2: Errors and convergence orders for (5.4i on the Dziuk surface with Choice 2 (analysis). 



Table shows the L? and DG errors for Choice 2. As expected, the experimental orders of conver- 
gence (EOCs) match up well with the theoretical convergence rates. Figure [T] shows the resulting DG 
approximation to (5.4i on the Dziuk surface using Choice 2. 

Figure 2(a,b) shows respectively the ratios of the L? and DG errors p^i with i ~ 1,3,4 where Err; 
denotes the error in the corresponding norm when using Choice /. Choices 2 (analysis) and 3 (average) 
appear to give the best results in both the L^ and DG norms. In particular, the additional symmetry 
induced by using Choice 3 which we mentioned previously makes it the preferable choice. 

A few remarks on Choice 4 with the true penalty term which, as mentioned before, would correspond 



to the Arnold et al. (2002 1 IP method on Ff,: interestingly, the scheme fails to converge for such a choice. 
The numerical scheme appears to be particularly sensitive to small perturbations in the off-diagonal 
entries of the resulting matrix, namely the ones caused by the product of the conormals n^ • n^ when 
using the true penalty term for Choice 4. Note that in the flat case, n,^ • nj^ is equal to — 1 . We tried 
to reproduce this problem in the flat case, taking two different values for the penalty parameter on e/, 
depending on whether we are assembling the diagonal or the off-diagonal block. Already a factor of 
10^^ leads to similar problems with stability. Since Choice 4 with or without the true penalty term was 
always less accurate than the other choices, we omit this choice in our next test case. 



5.4 Test Problem on Enzensberger-Stem Surface 

:{x e I 



Our next test problem considers (5.4 1 on F 
40 



: 400(xV+yV 



2 2 



,2^3 _ 



0} whose exact solution is again chosen to be given by m(x) — x\X2 
problem, there is no explicit projection map so we make use of the first order ad-hoc algorithm. In this 



)-(i-x2-r-z- 

As for the previous test 



test problem, the computation of Zij-M to derive the right-hand side of ( 5.4 1 is done via our approximation 



of the Laplace-Beltrami operator described in Section 5.1 
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(a) 



Fig. 1: DG approximation of (5.4 1 on the Dziuk surface with Choice 2 (analysis). 
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RaliD of L errors w r.t benchmark error for different normal choices (conforming) Ratio of DG errors w.r.t benchmari^ error for different normai cfiotces (conforitling) 
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(b) 



Fig. 2: Ratio of respectively L? and DG errors for (5.4i on the Dziuk surface with respect to 
the analysis error (Choice 2) for Choices 1, 3 and 4. 
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Elements 


h 


L2 -error 


L2-eoc 


DG-error 


DG-eoc 


2358 


0.163789 


0.476777 




0.998066 




9432 


0.0817973 


0.175293 


1.44 


0.472241 


1.08 


37728 


0.040885 


0.0160606 


3.45 


0.150144 


1.65 


150912 


0.020441 1 


0.00139698 


3.52 


0.0703901 


1.09 


603648 


0.0102204 


0.00033846 


2.04 


0.03473453 


1.02 


2414592 


0.00511 


7.86713e-05 


2.10 


0.0172348 


1.01 



Table 3: Errors and convergence orders for (5.4i on the Enzensberger-Stern surface with Choice 2 
(analysis). 



Table[3jshows the L^ and DG errors for Choice 2. Although the EOCs are more erratic than for the 
previous test problem, largely due to our approximation of the Laplace-Beltrami operator, they neverthe- 
less match up well with theoretical convergence rates. Figure l3] shows the resulting DG approximation 
to ( 5.4 1 on this surface using Choice 2. We again consider the DG approximation of ( 5.4 1 for different 
choices of n^, n+ and nj . Figure W(a,b) shows respectively the ratios of the L? and DG errors. These 
results confirm that Choices 2 and 3 are the preferable ones to use for DG schemes on surfaces. 



6. Extensions 

Although our analysis was restricted to conforming grids due to the nature of the surface approxima- 



tion, our numerical tests suggest that the estimates of Theorem 4. 1 also hold for nonconforming grids 
as shown in Table |4] for the Dziuk surface. Future work aims to derive a-priori error estimates for 



Elements 


h 


L2 -error 


L2-eoc 


DG-error 


DG-eoc 


230 


0.353599 


0.21889 




0.777436 




920 


0.176993 


0.0530078 


2.05 


0.413817 


0.91 


3680 


0.0885231 


0.0281113 


0.92 


0.223119 


0.89 


14720 


0.0442651 


0.00442299 


2.67 


0.111518 


1.00 


58880 


0.022133 


0.00104207 


2.08 


0.0562128 


0.99 


235520 


0.0110666 


0.00026444 


1.99 


0.0281247 


1.00 


942080 


0.00553329 


6.60383e-05 


2.00 


0.0140544 


1.00 



Table 4: Errors and convergence orders for (5.4i on the Dziuk surface with Choice 2 (analysis) for 
a nonconforming grid. 



nonconforming grids. 



Demlow ( 2009 1 has proven that in particular, for a linear approximation of the surface and quadratic 
polynomial basis functions, the FEM error scales quadratically in both the L? and //' norms. Numerical 
tests suggest that our DG scheme scales similarly in the L^ and DG norms as shown in Table [s] for 
the Dziuk surface. In future work we aim to derive higher order a-priori error estimates (that is, both 
higher order polynomial basis functions and higher order approximations of the surface) for the DG 
approximation in a similar fashion to the work done inlDemlow (2009[l. 
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(a) 



Fig. 3: DG approximation of (5.4 1 on the Enzensberger-Stern surface with Choice 2 (analysis). 
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(b) 



Fig. 4: Ratio of respectively L^ and DG errors for (5.4i on the Enzensberger-Stern surface with 
respect to the analysis error (Choice 2) for Choices 1 and 3. 
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Elements 


h 


L2-error 


L2-eoc 


DG-error 


DG-eoc 


92 


0.704521 


0.136442 




0.322416 




368 


0.353599 


0.0551454 


1.31 


0.150303 


1.10 


1472 


0.176993 


0.0215041 


1.36 


0.0601722 


1.32 


5888 


0.0885231 


0.00448861 


2.26 


0.0182412 


1.72 


23552 


0.0442651 


0.00120287 


1.90 


0.00513161 


1.83 


94208 


0.022133 


0.00029651 


2.02 


0.00130482 


1.98 


376832 


0.0110666 


7.41044e-05 


2.00 


0.00032728 


2.00 



Table 5: Errors and convergence orders for (5.4 1 on the Dziuk surface with Choice 3 (average) 
using quadratic polynomial basis functions. 
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